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ABSTRACT 

An innovative hyperbolic preconditioning technique is developed 
for the numerical solution of the Helmholtz equation which governs acous- 
tic propagation in ducts. Two pseudo-time parameters are used to pro- 
duce an explicit iterative finite difference scheme. This scheme eliminates 
the large matrix storage requirements normally associated with numeri- 
cal solutions to the Helmholtz equation. The solution procedure is very 
fast when compared to other transient and steady methods. Optimization 
and an error analysis of the preconditioning factors are present. For vali- 
dation, the method is applied to sound propagation in a 2D semi-infinite 
hard wall duct. 


INTRODUCTION 

The Helmholtz equation plays an important role in the study of 
acoustics as well as electromagnetic propagation and quantum mechan- 
ics. Unfortunately, large matrix storage requirements are generally asso- 
ciated with numerical solutions of the Helmholtz equation. To reduce 
these requirements. Bay less, Goldstein, and Turkel (1982) developed an 
iterative approach to solve the associated matrix equation. More recently, 
Baumeister and Kreider (1996) developed a preconditioned transient fi- 
nite difference scheme to solve the Helmholtz equation as well as the 
more general linearized potential flow equations. Their introduction of 
time dependence into the Helmholtz equation eliminates the large ma- 
trix storage requirements of the algorithm. This paper is concerned with 
the development of a more efficient preconditioning method to acceler- 
ate convergence for the Helmholtz equation. 

A standard technique of solving steady state partial differential equa- 
tions is to march their time dependent form until the steady state is reached. 
When the transient is not of interest, acceleration parameters can be 
employed to speed the convergence. This type of differential manipula- 
tion is often associated with preconditioning of both time dependent 
(Turkel, Fiterman and van Leer, 1993) and time independent (Turkel and 
Amone, 1993) partial differential equations. Generally, acceleration 
parameters destroy the time accuracy of the solution. Turkel (1982, pp. 31) 


addresses many of the issues associated with the acceleration to a steady 
state solution. 

In this paper, new preconditioning factors are introduced to speed 
the convergence of the transient finite difference scheme in solving the 
Helmholtz equation. Solution times are reduced by an order of magni- 
tude over the previous approach. For validation, the method is applied to 
plane wave sound propagation in a 2D semi-infinite hard wall duct. The 
paper contains a description of the problem, the brief development of the 
preconditioning technique, the introduction of the acceleration param- 
eters, the finite difference formulation with a stability analysis, and sev- 
eral numerical examples with error estimates. 


NOMENCLATURE 

C Q # dimensional speed of sound 

C dimensionless speed of sound, Eq. (2) 

D # dimensional duct height 

D duct height, D # /D # , D = 1 

e k total Lj convergence error at step k, Eq. (16) 

e^ e k error per axial wavelength, Eq. (17) 

f* dimensional frequency 

f dimensionless frequency, f*D # /C Q # , Eq. (2) 

i 

L length, L # /D # 

IMI absolute value of Mach number, Eq. ( 12) 
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n unit outward normal 
t dimensionless time, f*t # 


HELMHOLTZ EQUATION AND BOUNDARY CONDITION 

The governing differential equation for studying wave propagation 
can be formulated in terms of a potential as 


tj. dimensionless total calculation time 
At time step 


£,2 ^xx+^yy or f ^tt “ ^xx + 


yy 


( 1 ) 


x dimensionless axial coordinate, x # /D # 

Ax axial grid spacing 

y dimensionless transverse coordinate, y # /D # 


where <J>'(x,y»t) is the dimensionless potential and subscripts indicate par- 
tial differentiation with respect to subscripted variables. The conventional 
normalization factors used to develop these nondimensional equations 
are given in the NOMENCLATURE. The dimensionless frequency f is 
defined as 


Ay transverse grid spacing 
a acceleration parameter 


f = 


f # D# 


C 


(2) 


P acceleration parameter 

<t>' transient potential, <(> # /C 0 # D # Eq. ( 1 ) 

<(> transient potential in frequency domain, Eq. (6) 


where the superscript # indicates a dimensional quantity. 

There are several ways to develop a frequency domain formulation 
for Eq. (1). The Fourier Transform can be applied if the potential has a 
multi-frequency content. In the monochromatic case, this is equivalent 
to assuming that 


\|/ Fourier transformed potential, Eq. (3) 
co dimensionless frequency, 2?rf 


<|>'(x,y,t) = y(x,y)e ,t0#t# = y(x,y)e l2m (3) 

which transforms Eq. (1) to the Helmholtz equation 


Subscripts 

i axial index, see Fig. 1 
j transverse index, see Fig. 1 

o ambient or reference condition 

Superscript 

# dimensional quantity 


° = ¥xx+Vyy +“ 2 V (4) 

where co = 2 n f. 

At the entrance of the duct, x = 0, the source is assumed to have the 

form 

<(*' = e -i2, “ or (5) 

Also, the duct is assumed to be semi-infinite in length so that waves 
propagate only to the right. 


k time step 

complex conjugate 

PROBLEM STATEMENT 

The problem under consideration here is the development of pre- 
conditioning acceleration factors to obtain the solution to the Helmholtz 
equation. This method will have application to the general study of sound 
propagation in ducts. The goal of the paper is to develop a stable, explicit 
finite difference scheme that significantly reduces the computation time 
required to solve the Helmholtz equation with a monochromatic noise 
source. The formulation is applied to a semi-infinite duct with a planar 
source at the duct inlet, as shown in Fig. 1 . 


PRECONDITIONED HELMHOLTZ EQUATION 

The Helmholtz equation is preconditioned by assuming that 

4»'(x.y,t) = <Kx,y,t)e -I “ ‘ = <)>(x,y,t)e _l2m (6) 

See Baumeister and Kreider (1996) for further discussion. This differs 
from the classical monochromatic transformation in that the amplitude $ 
(no prime) is no longer independent of time. Under this transformation, 
Eq. (1) becomes 

f\ t - 2ifco<t> t = <|> xx + <|> yy + co 2 <J) (7) 
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The only difference between the Helmholtz Eqs. (4) and (7) is 
the presence of the time derivative terms on the left hand side. Physi- 
cally, the time dependence in <|)(x,y,t) is caused by assuming that the duct 
is quiescent at time 0, and that the source is turned on at that instant. A 
series of numerical calculations reported by Baumeister and Kreider 
(1996) show that 


lim <t>(x,y,t) = y(x,y) 

t— >oo 


when the f 2 <(> tt is dropped (parabolic approximation). 


( 8 ) 


ACCELERATION PARAMETERS 

To speed the convergence to the steady state solution \y(x,y), accel- 
eration parameters a and p are added to Eq. (7) as follows: 

-p2ifaxj) t = <() xx + <|>yy + co 2 c}> (9) 

This is a generalization of the preconditioning done in Baumeister and 
Kreider (1996), which dealt with the a = 0, p = 1 case (parabolic 
approximation). After the formulation of the difference equations, 
the acceleration of convergence is tested over a range of acceleration 
parameters. 


FINITE DIFFERENCE EQUATIONS 

The potential at the spatial grid points (xj,yp (Fig. 1 ) is determined 
by iterating the initial condition over time steps t k = kAt. Away from the 
duct boundaries, each partial derivative in Eq. (9) can be expressed using 
central differences, which yields 


k+lf af2 P itof 

:: Ut 2 At, 


iKI 

^i.j 


ik 

Ai 


f 2af 2 2 2_ 

At 2 Ax 2 Ay 2 




(10) 


. ik 

+ ^ij+l 


( i \ 


v A y' j 


+ ‘t'f.j-i 


( , 


v A y ) 


+ <j 


k-if af 2 picof 


At 2 At 


where Ax, Ay, and At are the space and time mesh spacings respectively, 
and = (Kx^yjA. Equation (10) is an explicit two step scheme. The 
field values at t° and t' 1 are assumed zero because the initial field is 
quiescent. 


The expressions for the difference equations at the hard wall bound- 
aries (y=0 & y = 1 ) employ the boundary condition 

V<))»n = 0 (11) 

where n is the unit outward normal. Baumeister (1980) gives precise 
details for generating the difference equations on the boundaries. 


STABILITY 

A von Neumann stability analysis is used to determine the condi- 
tions on At, Ax, and Ay required for conditional stability as a function of 
the acceleration parameters a and p. Conditional stability means that the 
amplification factor, which describes how errors propagate from one time 
step to the next, has magnitude one. Thus, when At, Ax, and Ay satisfy 
the stability criteria, errors are not magnified or diminished in magni- 
tude. This is a desirable property, since the numerical formulation can- 
not distinguish between a roundoff error and a small physical oscillation. 

For the case a = 0 and p = 1, treated by Baumeister and Kreider 
(1996) and herein denoted the parabolic preconditioner (because the sec- 
ond order time derivative does not appear), the stability analysis indi- 
cates that the method is conditionally stable, subject to the condition 
(conservative estimate of Mach number effects) 



In a typical application, f is set by the operating conditions in the duct. 
The grid spacing parameters Ax and Ay are set to resolve the estimated 
spatial harmonic variation of the potential field and At is chosen to sat- 
isfy Eq. (12). Of course, the stability analysis does not take into account 
boundary conditions. For stability, gradient boundary conditions gener- 
ally require the use of smaller At than predicted by Eq. (12). 

For the case a * 0 and P * 0, herein denoted the mixed preconditioner, 
the stability analysis indicates that the method is conditionally stable, 
subject to the two conditions 

At < -1 J(p 2 - a) (13) 


and 


At < 


af 

1 1 

- + - 


co 


Ax z Ay 


p7rf 2 


1 1 

- + - 


(14) 


to 


Ax z Ay 


Generally, the second squared term on the right side of (14) is much 
smaller than the first term. For a typical scenario (Ax = 0.05, Ay =0.1, 
f = 1, a = 0.95, p = 1), At < 0.071 from Eq. (13) and At < 0.044 from 
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Eq. (14). Consequently, the stability criteria yield At<.044 in this mixed 
case example. 

NUMERICAL EXAMPLES 

Let a plane wave propagate from the left into a semi infinite quies- 
cent duct. The potential field is to be computed in the duct for 0 < x < 1. 
Examples 1 to 3 illustrate the effects of changing a and p with dimen- 
sionless frequency f = 1 on the grid specified by Ax = 0.05 and Ay = 0.5. 
Example 4 shows the effect of increasing the frequency, by using f = 5. 

Because boundary conditions can introduce instabilities (Baumeister, 
1982 and Cabelli, 1982) into otherwise stable finite difference schemes, 
it is important to test the iteration scheme for convergence in the absence 
of an exit boundary condition. Therefore, in these examples, the compu- 
tational boundary is set at x = 50, far enough away from the true bound- 
ary x = 1 that any artifacts arising from the exit boundary condition do 
not affect the solution in 0 < x < 1 . 

The numerical results for plane wave propagation are compared to 
the exact results of the steady state solution, given by 


parabolic preconditioning because the number of iterations required for 
convergence is an order of magnitude lower. 

Example 3. Mixed Preconditioning — a = 16, p = 4.25, f = 1. 

The time increment At is set at 0.2. The maximum stable value, 
from Eq. (14), is 0.201. The error as a function of the number of itera- 
tions is shown in Fig. 5. After 500 iterations, e 100 = 0.0136225. The 
results here are nearly identical to those shown in Fig. 4. The numerical 
solution again shows excellent agreement with the analytic solution. As 
before, the plots of magnitude and phase are omitted because they are 
virtually identical to Figs. 2(a) and (b). In this example, with the intro- 
duction of the large acceleration parameters, the time variable loses its 
physical meaning. In effect, the large a reduces the effective speed of 
propagation C in Eq. (1) which then requires a longer time for the tran- 
sient <|> to approximate the steady state solution \| / as required by Eq. (8). 
However, the critical parameter is the number of iterations required to 
obtain a solution. The convergence rates for examples 2 and 3 are 
nearly identical. 


\|/(x) = e icox 


(15) 


Example 4. Mixed Preconditioning — a = 0.95, p = 1, f = 5. 


In the region 0 < x < 1, the Lj norm of the global error e k between the 
exact solution \\f and the numerical solution ^ at time step k is used as a 
meas-ure of the convergence. The error is defined as 

e k = J 0 J 0 J(«t )k - -v)dy dx. (16) 

Example 1. Parabolic Preconditioning — a = 0, P = 1, f = 1. 

The choice of a = 0 eliminates the second order time derivative. 
The time increment At is set at 0.007. The maximum stable value, from 
Eq. (12), is 0.00797. The numerical and exact solutions are compared in 
Fig. 2(a) (real and imaginary parts of the potential) and in Fig. 2(b) (mag- 
nitude of the potential) after 1000 iterations. The numerical solution shows 
excellent agreement with the analytic solution. The error as a function of 
the iteration number is shown in Fig. 3. After 1000 iterations, e 1000 = 
0.017 1 . This residual error is the result of the usual round off and trunca- 
tion errors associated with finite differences. 


To resolve the shorter wavelengths associated with the higher fre- 
quency f = 5, Ax is reduced by a factor of 5 to 0.01. The time increment 
At is set at 0.04. The maximum stable value, from Eq. (14), is 0.0497. 
The numerical and exact solutions are compared in Fig. 6(a) (real and 
imaginary parts of the potential) and in Fig. 6(b) (magnitude of the po- 
tential) after 1000 iterations. The numerical solution again shows excel- 
lent agreement with the analytic solution. The error as a function of 
iteration number is shown in Fig. 7. After 1000 iterations, e 1000 = 0.0659. 
This error is about 5 times higher than that for the f = 1 case. Since the 
two solution plots show roughly the same degree of accuracy, the change 
in total global error is apparently caused by the fact that there are 5 times 
as many grid points in the f = 5 case. The number of grid points is pro- 
portional to the number of axial wavelengths or frequency. So dividing 
by the frequency gives a rough measure of the local error at each grid 
point, which is, in this case, 

e kX=y = 0.0132 (17) 


Example 2. Mixed Preconditioning — a = 0.95, P = 1 , f = 1 . 


It should be noted that acceptable solutions can be obtained using fewer 
iterations than the asymptotic values indicated by Figs. 7 to 9. 


The time increment At is set at 0.04. The maximum stable value, 
from Eq. (14), is 0.049. The error as a function of the iteration number is 
shown in Fig. 4. After 100 iterations, the error is e 100 = 0.0136201. It is 
important to note that the number of iterations required for convergence 
drops by an order of magnitude with mixed preconditioning as com- 
pared to parabolic preconditioning (shown by dashed line in Fig. 4). 
Therefore, a calculation which took 1 minute with the parabolic approach 
would now take only 6 seconds with the Hyperbolic Approach. Again, 
the numerical solution shows excellent agreement with the analytic so- 
lution after 100 iterations. The plots of the magnitude and the real and 
imaginary parts of the potential are virtually identical to Figs. 2(a) and 
(b), and hence are omitted. Clearly, mixed preconditioning is superior to 


Optimal Acceleration Parameters 

Numerous numerical calculations were performed to determine the 
optimal choice of a and p to reduce the number of iterations to conver- 
gence. In these calculations a was set to 1 and p varied. Other choices 
of a yield approximately the same rate of convergence. For 
convergence with minimum iterations, the time increment At and p were 
set as follows: 


^opt 


af 

1 1 

- + - 


co 


Ax' 


Ay' 


(18) 
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( 22 ) 



a + rc 2 At 


2 

opt 


(19) 


t 


= — = Lf = 1*5 = 5 
C 


For a = 1, Ax = 0.05, Ay = 0.5 and f = 1 , the optimal time increment 

is calculated to be 1.01244 
the number of iterations is 
shown in Fig. 8. After 100 iterations, e 100 = 0.0137. The converged error 
is nearly identical to the a = 0.95 and (3 = 1 .0 nonoptimal case shown in 
Fig. 4. However, the optimal case reaches this error with nearly one order 
of magnitude in fewer iterations. Also, the error curve is seen to be much 
smoother. Therefore, a calculation which took 1 minute with the para- 
bolic approach would now take only 1.8 seconds with the Hyperbolic 
Approach. The numerical solution again shows excellent agreement with 
the analytic solution. As before, the plots of magnitude and phase are 
omitted because they are virtually identical to Figs. 2(a) and (b). 

The number of iterations to obtain the final solution can be reduced 
even further by increasing the spatial distance between nodes. For a = 1, 
Ax = 0.0833 (66% increase), Ay = 0.5 and f = 1, the optimal time incre- 
ment At t from (18) is 0.0851. The optimal P opt is calculated to be 
1.03614 from Eq. (19) with a 1.001 factor of safety. The error as a func- 
tion of the number of iterations is shown in Fig. 9. After 100 iterations, 
e l00 = 0.0398 which is approximately 3 times larger than the previous 
case shown in Fig. 8 with Ax = 0.05. However, this error is still accept- 
able for accurate solutions. In fact, after 16 iterations, the error level has 
dropped to 0.04732 which gives excellent results for phase and magni- 
tude plots of the potential as shown in Figs. 10(a) and (b). 

The speed at which the transient solution approaches the steady 
state Helmholtz solution is now discussed. 

Speed of Propagation 

For the conventional wave equation, Pearson (1953) has shown that 
the total time tj required for the steady state solution (15) to become 
established in the duct is equal to the time for a plane wave propagating 
at the speed of sound to reach the end of the duct. Similarly, with the 
preconditioned Helmholtz Eq. (9) and a = 1 , the steady state solution 
y(x) to the Helmholtz equation propagates outward at the dimensionless 
speed C, where from Eq. (2) 


At opt ,from (18), is 0.0503. The optimal (3 ( 


opt 


from Eq. (19). The error as a function of 



f #D# 


( 20 ) 


As seen in Fig. 12, the solution moves to the right with a distinct front at 
the speed C = 0.2. The front arrives at the exit at t ~ 5.0. 

The number of axial grid points and the number of iterations re- 
quired for convergence are both directly proportional to frequency. A 
comparison of Figs. 11 and 12 shows the factor of 5 difference in the 
number of required iterations when the frequency is increased by 5. For 
increased frequency, the solution moves more slowly to the right. Gener- 
ally, the solution time should be increased by 30 percent over that pre- 
dicted by Eqs. (21) or (22) for more accurate results. 


SOLUTION METHODS 

With the approach developed in this paper, three different solution 
techniques are now available to solve the Helmholtz equation, as shown 
in Fig. 13. The Fourier transform approach in the right column, with 
finite differences or finite elements, results in a matrix equation. Because 
this matrix is not positive definite, matrix elimination solutions are gen- 
erally employed, requiring extensive computer memory for high frequency 
propagation. The transient solution to the hyperbolic wave equation, 
shown in the left column, eliminates matrix storage requirements by 
iterating finite difference approximations to the steady state solution. 

The third option, preconditioning using the mixed formulation (a * 0, 
|3 * 0), is shown in the center column. This approach eliminates matrix 
storage requirements, and has less stringent stability conditions than the 
transient solution, so it converges in fewer iterations. 


CONCLUSION 

Accelerated numerical preconditioning of the Helmholtz equation 
has been developed. The field is iterated in time from an initial value of 
0 to attain the steady state solution. The method eliminates the large 
matrix storage requirements of steady state finite difference or finite el- 
ement techniques in the frequency domain. In each example provided, 
the numerical solution quickly and accurately converges to the exact 
steady state solution. The hyperbolic preconditioning developed in this 
paper has more than an order of magnitude faster convergence than the 
previously developed parabolic preconditioning approach. 


For a frequency f = 1, Fig. 1 1(a) shows the developing wave front 
(real and imaginary parts) as a function of the iteration number while 
Fig. 1 1(b) shows the magnitude of the potential. For this case, 

t = — = Lf = l*l = l (21) 

C 

As seen in Fig. 1 1, the solution moves to the right with a distinct front at 
the speed C = 1. The front arrives at the exit at t « 1 .0. 

Similarly, for a frequency f = 5, Fig. 12(a) shows the developing 
wave front (real and imaginary parts) as a function of the iteration num- 
ber while Fig. 12(b) shows the magnitude of the potential. For this case, 
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Figure 1 . — Structured finite difference-time dependent (FD-TD) 
mesh for semi-infinite rectangular duct. 
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Figure 2.— Analytical and numerical potential profile along wall for 
plane wave propagating in a semi-infinite hard walled duct (f = 1). 
(a) Real and imaginary parts, (b) Magnitude. 
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Number of iterations 


Figure 3. — Integrated solution error as a function of the number of 
iterations of the finite difference equations (f = 1) — ► parabolic 
approximation. 



Number of iterations 

Figure 4. — Integrated solution error as a function of the number of 
iterations of the finite difference equations (f = 1). 
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Figure 5. — integrated solution error as a function of the number of 
iterations of the finite difference equations (f = 1). 
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Figure 6.— -Analytical and numerical potential profile along wall for 
plane wave propagating in a semi-infinite hard walled duct (f = 5). 
(a) Real and imaginary parts, (b) Magnitude. 
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Number of iterations 


Figure 7. — Integrated solution error as a function of the number of 
iterations of the finite difference equations (f = 5). 
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Number of iterations 


Figure 8.— Integrated solution error as a function of the number of 
iterations of the finite difference equations (f = 1) with optimal 
parameters. 



Number of iterations 


Figure 9.— Integrated solution error as a function of the number of 
iterations of the finite difference equations (f = 1). 
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Figure 10. — Analytical and numerical potential profile along wall for 
plane wave propagating in a semi-infinite hard walled duct (f = 1). 
(a) Real and imaginary parts, (b) Magnitude. 
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Real and imag potential Real and imag potential Real and imag potential Real and imag potential 



Figure 1 1 . — Developing history of disturbance propagation in Fourier transformed domain as a 
function of number of iterations and time (f = 1). (a) Real and imaginary parts, (b) Magnitude. 
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Figure 12. — Developing history of disturbance propagation in Fourier transformed domain as a function 
of number of iterations and time (f = 5). (a) Real and imaginary parts, (b) Magnitude. 
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Figure 13. — Alternate finite difference/element methods in solving wave equations. 
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